load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

location_fc="/Users/ivana/WORK_BSC/data_to_publish/EC-Earth3P/control/"
location_fw="/Users/ivana/WORK_BSC/data_to_publish/EC-Earth3P/lowice/"

filsw = systemfunc ("ls "+location_fw+"toafluxes_allmons_ym_ensmean_y1130_tm.nc")
filsc = systemfunc ("ls "+location_fc+"toafluxes_allmons_ym_ensmean_y1130_tm.nc")

fw    = addfiles (filsw, "r")
fc    = addfiles (filsc, "r")

;print(fw)
;print(fc)

;*******************************************;
;Read in variables
;*******************************************
lat=fc[0]->latitude
lon=fc[0]->longitude
time=fc[0] ->time
rad  = 4.0 * atan(1.0) / 180.

clat=dble2flt(cos(lat(:,0)*rad))

ntime= dimsizes(time)

rsdtc=fc[:]->rsdt
rsdtw=fw[:]->rsdt

rsutc=fc[:]->rsut
rsutw=fw[:]->rsut

rlutc=fc[:]->rlut
rlutw=fw[:]->rlut

tempsc=rsdtc(0,:,:)
tempsc=rsdtc(0,:,:)-rsutc(0,:,:)-rlutc(0,:,:) 

tempsw=rsdtc(0,:,:)
tempsw=rsdtw(0,:,:)-rsutw(0,:,:)-rlutw(0,:,:) 

printVarSummary(tempsc)
printVarSummary(tempsw)


;tempsc!0 = "lat"
;tempsc!1 = "lon"
;tempsw!0 = "lat"
;tempsw!1 = "lon"

tempsc_z=tempsc(:,0)
tempsc_z=dim_avg(tempsc(j|:,i|:))
tempsw_z=tempsw(:,0)
tempsw_z=dim_avg(tempsw(j|:,i|:))
;************************
; calculate the anomaly:
;************************

  finfo_ano=tempsw_z
  finfo_ano=tempsw_z - tempsc_z

  printVarSummary(finfo_ano)
  printVarSummary(clat)

  finfo_ano=finfo_ano*clat

  finfo_ano!0="lat"
  finfo_ano&lat=lat(:,0)
  ;print(finfo_ano)

;________________
;write data

;name ="fluxnet_ano_zonal_ccsm4som"
;outf = addfile(name+".nc","c")
;outf->time = time
;outf->lat  = lat
;outf->finfo_ano   = finfo_ano

;****************************************
; Create plot
;****************************************
wks = gsn_open_wks("pdf","zonal_toanetflux_ecearth3p")

   plot    = new (1,"graphic")

   res                      = True          ; individual plot
   res@gsnDraw              = False
   res@gsnFrame             = False
   res@xyLineThicknessF     = 2.0
   res@xyMonoLineColor         =True
   res@xyLineColor         = "black"
   res@gsnYRefLine           = 0.0             ; create a reference line
   res@gsnAboveYRefLineColor = "red"              ; above ref line fill red
   res@gsnBelowYRefLineColor = "blue"             ; below ref line fill blue

   res@xyDashPattern        = 0                  ; Make curves all solid

   res@vpHeightF= 0.4                    ; change aspect ratio of plot
   res@vpWidthF = 0.8

   res@trYMinF=-1
   res@trYMaxF=1.5

   res@tiXAxisFontHeightF = 0.010
   res@tiYAxisFontHeightF = 0.015

   res@gsnLeftString = ""
   res@tiYAxisString = "TOA net flux ano"
   res@tiXAxisString = ""

   plot(0)  = gsn_csm_xy (wks,lat(:,0),finfo_ano,res)

;***** make panel *****


   resP                     = True          ; panel resources
   resP@gsnMaximize         = True

   resP@txFontHeightF = 0.015
   gsn_panel(wks,plot,(/1,1/),resP)

